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STABILITY OF MULTISTEP METHODS IN NUMERICAL INTEGRATION 


By Robert N. Lea 
Manned Spacecraft Center 


SUMMARY 


One widely used technique for the numerical solution of a differential 
equation is to approximate the differential equation by a difference equation, 
and then to solve this difference equation. This paper treats the stability 
of the solutions obtained by this method. An original development leading to 
a definition of stability shows a relation between stability and certain prop- 
erties of the systems of differential equations to be solved. Previous investi- 
gations (with the exception of Dahlquist r s work) have shown a similar relation, 
but were limited to the consideration of a single equation. Dahlquist's theory, 
while valid for systems of equations, shows no such relation. 


INTRODUCTION 


A large part of the computer time at NASA Manned Spacecraft Center, 
Houston, Texas, is devoted to the numerical solution of differential equations. 
In order to reduce this computing time, it is necessary that accurate, more 
rapid methods of numerical integration be developed. One important criterion 
of accuracy in selecting a method is that it be stable. The stability of a 
method is that property which causes an error introduced at some step to tend 
to decay in succeeding steps, rather than to grow and eventually destroy the 
usefulness of the approximate solution. It is, therefore, extremely important 
that users of numerical integration procedures be familiar with the concept of 
stability. The purpose of this paper is to discuss stability of numerical 
methods, and to extend the definition of stability to systems of equations. 
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STABILITY AM) CONVERGENCE 


Consider the difference equation 

\ Y rH-k + VlU-l + • • • + Vn = h ( l 

to solve the ordinary differential equation 


b F + b n F . , + ...+ b.F (1) 

k n+k k-1 n+k-1 0 nj v 


Y' = F(x, Y) 

X(a) = Y q 

Let Y(x) be the solution to equation (2), and define the operator 

l[y(x)J = £ JsuYCx + ih) - hb^'Cx + ih)J 

By expanding L in a power series about x, it can be shown that 


( 2 ) 


L [y(x)J = O^ 1 ) 


•2 



if, and only if, the following p+1 linear relations hold: 


Statement (a) 


Statement (b) 


Statement ( c ) 



i=0 

k 


& 


a . = 0 
1 



b i fSrr] 



= 0 


(s — 2 , ••• , p) 

(s = 1) 


Definition A. Let p be the largest value of s for which 
statements (a) to (c) of the above assertion hold. Then p will be called 
the degree of the operator L or the degree of the difference equation (l). 


Definition B. The number of preceding points occurring, explicitly or 
implicitly, in equation (l) will be called the order of the operator L or 
the order of the difference equation (l). 

(The above definitions are those used by Dahlquist (ref. l). Some 
authors define the quantity of definition A as the order of the method and 
refer to the quantity in definition B as the step number of the method. ) 


As a simple example, consider the point-slope difference equation 


n+1 


- Y 


n 


hF 

n 


This method has order 1 since it uses only one preceding point to 

compute Y The equation 

lTy(x)J = Y(x + h) - Y(x) - hY'(x) 

is the operator associated with the difference equation. Statement (a) of defi- 
nition A holds, since a 1 = 1 and a Q = -1. The largest value of s for which 

statement (b) is true is determined in the following manner: 


If s = 1, then 


E 

i=0 


a . i 
i 


- b. 

i 


- 1 + 1=0 


and, if s = 2, 


then 


f a i .2 

E T 1 

1=0 


b. 1 
x 



Therefore, the point-slope formula has degree 1. 
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The polynominals 


P(?) = + + ••• + a Q 


cr( I ) = \ + Vi ^" 1 + " * + b 0 


(3) 


can be associated with the difference equation (l), and hence , with every 
pair of polynomials of the type of equation ( 3 ), can be associated a difference 
equation of the form of equation (l), provided the degree of p is at least 
equal to the degree of a. The following assumptions can be made concerning 
the polynomials of equation ( 3 ): 


Assumption (a). 
Assumption (b). 
Assumption (c). 


The coefficients a.. b.(i = 0, 

The polynomials p(|) and a( £ ) 

ihe degree p of the operator 
(Condition of consistency) 


... , k) are real, / 0. 
have no common factor. 

L is at least 1. 


Definition C. The difference equation (l) is said to be stable if the 
following assumptions are satisfied in addition to assumptions (a) to (c): 


Assumption (d). 
circle. 

Assumption (e). 


The roots of p(|) are located within or on the unit 
The roots on the unit circle are distinct. 


Theorem I. A necessary and sufficient condition for convergence of the 
linear multistep method of equation (l) is that the condition of stability is 
satisfied. 


The proof of Theorem I and also the proofs of the following three theo- 
rems are given in references 1 and 2. 

Theorem II. The degree p of a stabxe' operator of order K can never 
exceed K + 2. If K is odd, the degree cannot exceed K + 1. 

Theorem III. If an operator of even order K is stable, then the 
conditions 


and 



00 


are necessary and sufficient in order that it should be of maximum degree 
K + 2. All roots of p(|) then have unit modulus. 


Definition D. The formula of equation (l) is said to be open if b^ = 0, 
k 


and closed if b n _ f 0. 



Theorem IV. If the degree p is greater than the order K and the 
method is stable, then the difference equation must be closed. 

Actually b^/a^ * 0* It is possible to construct an open, stable operator 

of degree p = K, starting from any polynomial p(|), and satisfying 
assumptions (a) to (e). 

As an example of an unstable method, consider the difference equation 

Y ~ 2Y + Y n = F (h,x,Y) ( 5 ) 

n+1 n n-1 l v 7 ’ v 

to solve the initial value problem Y* = 0, Y(0) =0. 

The associated polynomial p(£) has a multiple root at | = 1 and 
therefore violates the condition of stability. Assume that Y ^ = Y ^ = 0 

are given exactly. Further, assume that an error is made in computing Y^ 

and the result is Y Q = c / 0. Using the formula of equation ( 5 ) and computing 

several values give the following results: 

Y 1 = 2e 
Y 2 = 3e 

Y = he 

3 

Y = (n + l)e 
n 

It is then seen that in an unstable method, any error introduced can 
continue to grow in succeeding calculations until the approximate solution 
becomes worthless. A stable method does not allow growth of this type. For 
example, consider the point-slope formula 

Y t - Y = hF (6) 

n+1 n n v 

and the more complicated formula 

Y n + 1 - k - k-1 - (7) 

The associated polynomials of equations (6) and ( 7 ) are respectively, 

p - 1 = 0 (8) 

and p 2 - -|p - | = 0 (9) 
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The root of equation (8) is p = 1 and the roots of equation ( 9 ) are 

p = 1 and p = Using equation (6) to solve the same problem as in the 

preceding example, it can be seen that the error in each step, after the error 
was introduced in the computation of Y^, has the same magnitude as the error 

in Yq. Using equation ( 7 ), a few values are computed: 


e 


2 


_ 



Notice that, assuming no further error is made, the error in step n is 
less than the error made in computing Yq. 

As Theorem I states, the condition of definition C is necessary and 
sufficient for convergence. 


WEAK INSTABILITY 


The second type of instability may allow unfavorable error growth in the 
immediate vicinity of zero step-size. Methods in which this can happen are 
called weakly or conditionally stable. This section is concerned with the 
integration of the single first-order equation 

y f = f(x,y) (10) 


using the difference equation (l). 

Define e. by the equation 

K 

f k ' r(S) - y k 

and e ’ by the equation 
K. 

s k ■ f [v y (*k)] - f k 
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It can 

be 

V ^ 

by the mean 


is 


be assumed that 
constant in the 
value theorem 


h is sufficiently small that, given a point 
interval with end points y^x^ and y^. Then, 


f df 

€ k = ^ 6 k 


(id 


By substituting the tine values into equation (l), and then the computed 
values, and subtracting the two, the following difference equation is obtained 
for the error e which was committed because the computed yalues were used 
rather than the true values. 


Vn+k + Vl'n+k-l + ••• + a 0 e n * h (Vwk + b k-l e Mk-l + ••• + VA) 

( 12 ) 

^f 

On using equation (ll) and defining h = h, equation (12) can be written as 


(=k - hb k) e n+k + (\-l - hb k-l) Vk-1 + ' ' ' + ( a 0 ' hb o) e n ' 


Now 


6 n = C l pn + C 2 P 2 n + ' * ‘ + 


(13) 


(1^) 


where the p^,i =!>••• > k are the roots of the associated polynominal 


(\ - hb k) pk + (vi - h \-i) plt ' 1 + • • • + C a o - hb o) = 0 (15) 

It should be noted that if the roots of equation (15) are not simple, the form 
of the previous equation for (eq. (l4)) will be slightly different. As an 

example, if is a root of multiplicity 2, 


E n = ( C 1 + C 2 n ) P i n + 


n , n 

C 3 P 3 + • ' • + Vk 


Relative stability can now be defined as follows: 


Definition E. In the region h < = 0, the difference equation (l) is said 
to be stable if the roots of the polynomial of equation (15) are contained 

c3f 

within the unit circle. Note that when h is positive, or is positive, 

the error already grows exponentially from the nature of the error equa- 
tion (ll). This growth is not regarded as instability in the method. 

The following theorem from reference 3 gives a stability criterion for 
a difference equation. 
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and 

A r,s \>s-l, k-r-1 

Thus , the number of computations involved in the application of Theorem V to 
a difference equation is reduced considerably. 

References 5 and 6 use the fact that the roots of a polynomial are con- 
tinuous functions of the coefficients of the polynomial to study stability of 
integration formulas. This property of roots of polynomials is proven in 

reference 7- It is possible to vary h in equation (15) and plot the roots 

as a function of h to determine the region of stability of equation (l). 

It is thought that the method of Theorem V is more precise than the afore- 
mention method; however , computing the characteristic roots of A might 

r > s 

not be simple. The previous method, while not quite as precise, should give 
a satisfactory bound on h. 

It is shown in reference 8 that if a predictor-corrector method is used 
where the corrector formula is used only once per integration step, the sta- 
bility of the method is dependent on the predictor formula as well as on the 
corrector. Therefore, in arriving at equation (15), an error equation must 
be derived for the predicted value, and this value must be substituted into 
the error equation for the corrected value. 

Finally, it should be noted that this latter concept of stability has not 
been extended to a system of differential equations. 
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EXTENSION OF THE DEFINITION OF STABILITY 


= Y^x^ - Yj and Z' = pjjxj , Y^x^j - F^. Substituting Y^x^ 


Let Z. = 

J 

and Y. into equation (l), and subtracting 'in the same manner as in the pre- 
J 

ceding section , yield the following difference equation for the error Z, 


a k Z n+k + 


+ a 0 Z n h k Z n+k 


+ V 


i) 


(16) 


Let 


where 


ij 


df . 

1 

^7 


A = a 


, the terms f. and y. "being the ith and jth components 


of the vectors P and Y, respectively. The assumption is now made that 

a. . is constant and it follows immediately from the mean value theorem that 
J 


z: = az. 

3 3 


( 17 ) 


Upon substituting equation (17) into equation (l 6) 9 the following differ- 
ence equation is obtained 


el Z ,+...+ a^Z = hA(b Z _ + . . . + b n Z \ 

^ n+k On \ k n+k 0 n ) 

For purposes of simplif ication, the following theorem is used. 

Theorem VI. Every complex nxn matrix A is similar to a matrix 
of the form 


/ 


j = 


0 

0 


0 \ 

0 


(is) 


(19) 


\ 0 


J S / 
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where 


and 


Jq is a diagonal matrix with diagonal X^, 

0 0 ... 




j. = 


X ,. 1 

q+i 

0 




Vi 


o 

0 


1 0 

• • • 

0 0 

0 0 


• • 


\ 


q+i 


1 

\+i/ 

(1=1, 


• , s) 


The 
need not 


\ > j — 1 j 2 y ••• 
J 

all be distinct. 


, q + s are the characteristic roots of A and 

If the X. are all simple, then A is similar to 
J 


a diagonal matrix whose diagonal elements are the characteristic roots of A. 
(For a proof, reference could be made to ref. 9 or almost any of the numerous 
texts on linear algebra and matrix theory. ) The form of equation (19 ) is 
noted by j(A). 


By the preceding theorem, there exists a non-singular matrix P, such 

that PAP" 1 = J(A). Set e = PZ, solve for Z = P" 1 ^, and substitute into 
equation (l8) to obtain 


a. e + 
k n+k 


a n € 

Q n 


= h 


'w(v.* 


V 


n) 


( 20 ) 


Transforming equation ( 17 ) in a similar manner, the following equation is 
obtained: 


«’ = J(A)e 


( 21 ) 


Only those error components of equation (21) which involve characteristic roots 
of A with negative real parts are considered, since the error in the other 
components already grow exponentially. This error growth due to the nature of 
the error equation (21) is not regarded as instability of the method. 


Let 


a, e . , , + 

k j , n+k 


v J; 


n 


- h \i( 


^k e j,n+k + 


b_e . \ 

0 J,nj 


+ hcfb, e . L , , , + 
V k j+l,n+k 


+ b 0 e J+l,n) (22) 
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be an arbitrary equation from the system of equations (20) where X . has 

J 

negative real parts. Note that c is either 0 or 1. If c is 0, only the 
roots of the following characteristic polynomial 

K - hx A) * k + • • • + (» 0 - h Vo) ’ 0 (23) 

need to be proved to lie within the unit circle. 


For the case c = 1, the following theorem will be needed. 

Lemma I. Let f(E) = a^ + a^E + . . . + a^E^ he a P^ynoniial with 
root R / 0. There exists an integer j ^ k such that 


£ a.i^R 1 ji 0 
i=0 


Proof. 


Consider the K + 1 homogeneous equations in k + 1 unknowns 

2 k ^ 

x q + Rx^ + Rxg + ... + R = 0 


Rx + 2R Xg + 


. . + kR 




= 0 


Rx^ + 2 k R^x^ + • . . + k^R^x^ = 0 


( 24 ) 


J 


The determinant of the matrix of coefficients of equation (24) reduces to 
the determinant 


R 

2R 2 . . 

. kR k 

R 

2 2 R 2 . . 

. k 2 R k 

R 

• • • • m 

2 k R 2 . . 

. kV 


11 



Factoring 
property that 


iR 1 



out of the ith column of equation ( 25 ) and using the 
|a T | for any square matrix A, equation (25 ) becomes 


1 1 
1 2 

1 k 



k ^ 1 


= c(-D k yy (i - j) 1 0 

l^i^^k 


(26) 


Since the determinant of the matrix of coefficients is non-zero, 
equation (24) can have only the trivial solution and therefore the coefficients 
of f(E) cannot be a solution. 


Theorem VII. Let E be an operator defined by 


Ey n - 


n+1 


and let 


f(E)y n = An r R n 


(27) 


be a non-homogeneous difference equation where R / 0 is a root of 

f(E) = a^ + a^E + . . . + a^E^ . Then, there exists an integer j, and 

r + 1 numbers, B. , B. B.., so that 
9 J J+l j+r' 

^ = Z B j + i nJ+iRn W 

i=0 


is a particular solution to equation ( 27 ). 


Proof. By Lemma I, there exists an integer j, so that 

k . . 

£ a^i^R 1 t 0 

i=0 

p 

Assume that j is the minimum such integer and that y^ has the form of 
equation (28). 
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Then 


* • | M 

f(E)y n P = R n £ aJ * 1 |B.(n + i )' 3 + B J+1 (n + i ) J+1 + • • • + B j +r ( n + l3J+r J 

i=0 (29) 


Equating coefficients of n 1 , i = 0 , . . . , r , yields the following 

r + 1 linear equations in r + 1 unknowns: 




k i=0 


k 

E 

i=0 


E a i lJE ‘l b j + I E a i iJ+lRl | Vi + -• + E 1- ~i 


a.i^Vl B. x =0 




i=0 


fj+l\ k 4 i 

E a i i<]R 

j / i=0 


B , + ... + 
0+1 


'j+r 


E a x lJ+r " lRi 


j+r-1 / i=0 


B = 0 
j+r 


) ( 30 ) 


[j+r\ k ± 

E "i 1 E 

J / i-0 


B. = A , 
J+r 7 


The coefficients of n 1 for i > r are identically zero. For considering 


r+m 


the coefficient of n , l^m^j , 


r+m+lA 


y a.R B + 

/ ) 1 r+m 

i=0 \ 1 


E a i iRB r+nH-l + ••• + 
i=0 


f r+j \ k . 

E 

V j-m / i =0 


Each B has a coefficient with a factor ^ a 1 i m R 1 where m<j. 

i=0 

Therefore the coefficient of each B is zero. 

Now consider again the non -homogeneous system of equations (30). It 
obviously has a solution, since its matrix of coefficients is triangular with 
all non-zero diagonal elements. Therefore , the r + 1 numbers 
B.,B. n . . . . , B. exist and the Theorem is proved. 

J+l J+r 

Let E be the operator defined in Theorem VII. That is. 


Ee i,n " € i,n+l 


13 



and 


f (E) = (a* - hX^E 11 + ( Vl - hX J b k _ 1 ) E k - 1 + ... + (a Q - bXjb 0 ) (51 ) 

where the a and b , 1=1,. . . , k, are defined by equation (l). 

Consider now equation (22) where c = 1. There exists some positive 
integer r such that 

A, . , — A* . 'i = • • • A# . 

j+r J+r-1 j 


and 


f (E)e = 0 

J+r,n 


(32) 


Assume further that r is the least such integer. Therefore, 


£ • , 
J +r 


m °s 

> A l-l n 
s A . n X 

Lj 1 S 

s=l i=l 


(33) 


where x. is a root of equation (23) of multiplicity j., i = 1, . . . , m. 


Substituting equation (33) into the equation for the following 

non-homogene ous difference equation is obtained: 


m s 

' Z E V x s" 

s=l i=l 


Therefore, for some integer p 


m p 

6 j+r-l,n 1 

s=l i=0 


1 

c .n x 
1 s 


n 


(3M 


(35) 


since the particular solution is of this form by Theorem VII and the solution 
of the homogeneous equation has the same form as equation (33)* Hence , by 
induction , there exists an integer q such that 


m 

S J 5 n " l 

S=1 


y D.ntc n 

L 1 s 

i=0 


(36) 
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Therefore, if x <1, s - 1, 


m. 


the error vector e in the transformed 


equations decreases to 0 as n increases and hence the error vector Z also 
converges to 0. Stability is now defined as follows: 


Definition F. 


Sf . 


Let X . , j =1, 
J 


, m, denote the characteristic roots of 


A = 


( a ij) 


, a^ . = The difference equation (l) is stable if, for all X^ 


iJ 


with negative real parts, the roots of 


(S - “V*) x + 


, a o ' *Vo 


)- 


= 0 


( 37 ) 


lie within the unit circle. 


CONCLUSIONS 


It has been shown that the region of stability of a method, if one exists, 
is dependent on the characteristic roots of the matrix A. If this information 
is available, it is possible to determine the step-size to insure that calcula- 
tions are within this stability region. If, however, there is interest only in 
the existence of a stability region, this can be guaranteed by requiring that 
all the roots of the associated polynomial, with the exception of the principal 
root, lie within the unit circle. 

Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, March 1, 1965 
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